import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt

# 设置中文字体，解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题


def extract_time_domain_features(signals):
    """
    提取时间序列信号的时域特征

    参数:
    signal: 一维数组，输入的时间序列信号

    返回:
    features_dict: 字典，包含提取的各个时域特征值
    """
    # 确保输入为numpy数组
    signals = np.array(signals)

    # 有量纲统计特征
    max_value = np.max(signals)  # 最大值
    min_value = np.min(signals)  # 最小值
    peak_to_peak = max_value - min_value  # 峰峰值
    mean_value = np.mean(signals)  # 均值
    abs_mean_value = np.mean(np.abs(signals))  # 绝对平均值
    std_value = np.std(signals, ddof=1)  # 标准差
    var_value = np.var(signals, ddof=1)  # 方差
    rms_value = np.sqrt(np.mean(signals ** 2))  # 均方根值(RMS)

    # 方根幅值 (需注意处理可能存在的负值)
    signal_abs = np.abs(signals)
    square_root_amplitude = (np.mean(np.sqrt(signal_abs))) ** 2

    # 高阶统计量
    skewness = stats.skew(signals)  # 偏度
    kurtosis = stats.kurtosis(signals)  # 峭度

    # 无量纲特征
    crest_factor = max_value / rms_value if rms_value != 0 else 0  # 峰值因子
    impulse_factor = max_value / abs_mean_value if abs_mean_value != 0 else 0  # 脉冲因子
    shape_factor = rms_value / abs_mean_value if abs_mean_value != 0 else 0  # 波形因子
    clearance_factor = max_value / square_root_amplitude if square_root_amplitude != 0 else 0  # 裕度因子

    # 将特征存入字典
    features_dict = {
        'max_value': max_value,
        'min_value': min_value,
        'peak_to_peak': peak_to_peak,
        'mean_value': mean_value,
        'abs_mean_value': abs_mean_value,
        'std_value': std_value,
        'var_value': var_value,
        'rms_value': rms_value,
        'square_root_amplitude': square_root_amplitude,
        'skewness': skewness,
        'kurtosis': kurtosis,
        'crest_factor': crest_factor,
        'impulse_factor': impulse_factor,
        'shape_factor': shape_factor,
        'clearance_factor': clearance_factor
    }

    return features_dict


# 示例用法
if __name__ == "__main__":
    # 生成示例信号（含噪声的正弦波）
    t = np.linspace(0, 1, 1000)
    signal = 0.5 * np.sin(2 * np.pi * 5 * t) + 0.1 * np.random.normal(size=len(t))

    # 提取特征
    features = extract_time_domain_features(signal)

    # 打印特征值
    print("时域特征提取结果:")
    for key, value in features.items():
        print(f"{key}: {value:.4f}")

    # 可视化信号和部分特征
    plt.figure(figsize=(12, 8))

    # 绘制原始信号
    plt.subplot(2, 1, 1)
    plt.plot(t, signal, label='原始信号')
    plt.axhline(y=features['mean_value'], color='r', linestyle='--', label=f"均值: {features['mean_value']:.3f}")
    plt.axhline(y=features['mean_value'] + features['std_value'], color='g', linestyle=':', label=f"均值±标准差")
    plt.axhline(y=features['mean_value'] - features['std_value'], color='g', linestyle=':')
    plt.fill_between(t, features['mean_value'] - features['std_value'],
                     features['mean_value'] + features['std_value'], alpha=0.1)
    plt.title('时间序列信号')
    plt.xlabel('时间')
    plt.ylabel('幅值')
    plt.legend()
    plt.grid(True)

    # 绘制特征值条形图（无量纲特征）
    plt.subplot(2, 1, 2)
    dimensionless_features = {k: v for k, v in features.items()
                              if k in ['crest_factor', 'impulse_factor', 'shape_factor', 'clearance_factor', 'skewness',
                                       'kurtosis']}
    plt.bar(range(len(dimensionless_features)), list(dimensionless_features.values()))
    plt.xticks(range(len(dimensionless_features)), list(dimensionless_features.keys()), rotation=45)
    plt.title('无量纲时域特征')
    plt.ylabel('特征值')
    plt.grid(True, axis='y')

    plt.tight_layout()
    plt.show()